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Many resources, such as oil, gas, or water, are extracted from porous soils and their exploration is often shared 
among different companies or nations. We show that the effective shares can be obtained by invading the porous 
medium simultaneously with various fluids. Partitioning a volume in two parts requires one division surface 
while the simultaneous boundary between three parts consists of lines. We identify and characterize these 
lines, showing that they form a fractal set consisting of a single thread spanning the medium and a surrounding 
cloud of loops. While the spanning thread has fractal dimension 1.55 ±0.03, the set of all lines has dimension 
1.69 ± 0.02. The size distribution of the loops follows a power law and the evolution of the set of lines exhibits 
a tricritical point described by a crossover with a negative dimension at criticality. 



I. INTRODUCTION 



II. MODEL 



Space partitioning is of interest in a wide spectrum of fields, 
ranging from materials science to medicine, with special rel- 
evance to computer graphics and the exploration of natural 
resources stored in soils. For example, if different compa- 
nies want to explore an oil reservoir they are interested in 
determining the volumetric share corresponding to each one 
inside the ground [1]. An additional degree of complexity 
comes into play when water is injected into the soil to push 
the oil to enhance extraction [2, 3]. Also in medical imag- 
ing, three-dimensional computed tomography scans need to 
be segmented to identify the different tissues. These pictures 
are discretized into pixels and a number is assigned to the 
bond between neighboring pixels corresponding to the inten- 
sity gradient. The resulting structure is similar to the one of 
a porous soil. By aggregating pixels pairwise from the low- 
est to the highest gradient it becomes possible to identify the 
boundaries of tissues [4]. 

Both problems consist in dividing space into parts: either 
the shares of the companies in the oil field or the different tis- 
sues in the image processing. In both cases, regions are sep- 
arated by division surfaces. Here we consider three regions 
and find that their division surfaces join in a fractal thread 
that crosses the medium, being surrounded by a cloud of dis- 
connected loops (see Fig. 1). In the case of oil exploration 
these points, where all three division surfaces merge, are the 
places where water should be injected to assure that no oil is 
pushed out on the wrong side. In medical image processing, 
the simultaneous boundary between three parts might indi- 
cate, for example, the region where a tumor is nested between 
two other tissues. 



A. Random media 



We consider a random medium consisting of pores arranged 
in a simple-cubic lattice connected through channels. To each 
channel k a threshold is randomly assigned following a uni- 
form distribution in the interval [0, 1). The fraction of open 
channels is tuned by a parameter p, such that channels with 
Pk < P are open while all the others are closed. Hereafter we 
use the language of fluids where p would correspond to the 
fluid pressure. For digital images, the pores would correspond 
to the pixels and the thresholds p^ to the intensity gradient 
between pairs of neighboring pixels. 
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FIG. 1 . Set of lines on which all division surfaces between three 
parts are in contact for a typical random medium. In addition to the 
backbone spanning the medium from left to right (shown in red), the 
set also contains a cloud of disconnected loops (green). The trans- 
parent planes are guides to the eye. 
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FIG. 2. Illustration of the model, (a) In the initial state (p = 0) the 
vertical faces of the cubic lattice are divided into three sets, (b) Ex- 
ample of the final state of the invasion (p= 1), dividing the medium 
into three parts: R, G, and B. RGB edges and nodes are shown as 
thick black lines and spheres, respectively. In (c) we separate the 
three pieces to be able to look inside. Solid lines represent the edges 
of the dual lattice of the pore network. The color of each cube cor- 
responds to the one of the fluid in the pore at the center of the cube. 
The channels connecting the pores are perpendicular to the faces of 
the cubes and for clarity they are not shown. The RGB edges and 
nodes are part of the dual lattice. 



B. Partitioning 

To find the partitioning of the medium into three parts, we 
divide the (four) vertical boundaries of a cubic system in three 
parts of about the same area. Each part corresponds to a differ- 
ent invading fluid distinguished by dyeing them with different 
colors: red (R), green (G), and blue (B) [see Fig. 2(a)]. In 
the illustration of Fig. 2 we see a medium of 5 x 5 x 5 pores. 
The pores are in the center of each cube and the edges are the 
bonds of the dual lattice of the pore network. The cubes have 
the color of the fluid contained in the corresponding pore. We 
invade the system simultaneously from all vertical walls. 

Starting with p = (i.e., all channels closed), the channel 
with the lowest threshold, in the invasion front, is selected and 
the fluid pressure p is increased to the value of this thresh- 
old. This channel and the empty pore connected to it are then 
invaded and colored according to the type of fluid that pene- 
trated into it. After that, invasion also cascades into all pores 
connected to this pore through channels with thresholds lower 
than the actual p. This process is repeated until all pores are 
invaded, under the constraint that the fluids can not displace 
each other, which does not allow to invade any pore by more 
than one fluid. 

In the final state, the medium is divided into three parts cor- 
responding to the pores filled either with an R, G, or B fluid. 
These parts are the maximum oil shares that each company 
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FIG. 3. Mass scaling analysis. The total mass M tot (□, total number 
of RGB nodes) and the mass of the spanning cluster M sc (0» number 
of nodes in largest RGB cluster) are shown as function of the lattice 
length L. One observes that for large lattices the masses scale as 
power laws of the lattice size, i.e., M to t ~ L dm and M sc ~ L dsc . The 
fractal dimensions are d tot = 1.69 ±0.02 and d sc = 1.55 ±0.03. In 
the inset we see the area of the division surfaces M sur (A) as function 
of the lattice size L. The fractal dimension of the division surfaces is 
dsw = 2.49 ± 0.02. Straight lines are guides to the eye. 



could extract from the exploration regions. This division is 
solely determined by the distribution of local thresholds, thus 
being intrinsic to the medium (for algorithmic details see Ap- 
pendix A). Here we will mainly focus on the final partitioning 
of the medium, corresponding to p = 1 (all pores invaded). 
However, we later also discuss the pressure p t at which two 
different fluids start to form an interface. 



C. RGB set 

An example for the partitioning into three parts is shown in 
Fig. 2(b). To better visualize the partitioning, we separate the 
three parts in Fig. 2(c). Every face that separates two colors 
is part of a division surface. Edges are shared by four differ- 
ent cubes. If three of the cubes sharing a common edge have 
different colors, we call this an RGB edge (thick black lines 
in Fig. 2). In fact, all RGB edges are on lines where all three 
surfaces separating regions of different color meet. The ver- 
tices attached to RGB edges are the RGB nodes and we define 
every set of nodes connected through RGB edges as an RGB 
cluster. All RGB nodes and edges of a medium form its RGB 
set. 



III. RESULTS 

The surfaces dividing two colors are fractal objects of frac- 
tal dimension d sur = 2.49 ± 0.02, as seen in the inset of Fig. 3, 
numerically consistent with what has been reported for wa- 
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tersheds in three dimensions [5, 6]. While these boundaries 
are singly connected, the RGB set consists of one spanning 
cluster connecting the two sides of the system surrounded 
by a cloud of smaller disconnected loops (see Fig. 1). As 
shown in Fig. 3, the entire RGB set is fractal of dimension 
dtot = 1.69 ±0.02, while the spanning cluster has a smaller 
fractal dimension d sc = 1.55 ±0.03. To analyze the topology 
of the spanning cluster, we used the burning method proposed 
in Ref. [7]. We found that the spanning cluster has loops, how- 
ever its backbone, elastic backbone, shortest path, and its set 
of singly connected RGB edges all have fractal dimensions 
consistent with d sc . 

The difference between d tot and d sc is due to the cloud of 
disconnected loops. These loops result from the entanglement 
of three compact regions, as illustrated in Fig. 4, which shows 
a transversal cross section of a medium, where the three re- 
gions are simultaneously in contact at different locations. In 
this particular case, the lower location (dashed circle) is where 
the spanning cluster intersects the cross section. The upper 
location (dotted circle) shows the cut through a disconnected 
loop: The G and B regions are in contact in an area completely 
surrounded by the R region, thus the contact line between the 
three forms a closed loop (discretization effects are discussed 
in the Appendix). The size distribution of the loops is shown 
in Fig. 5(a), where the size s is defined as the number of RGB 
nodes forming the loop. A power-law distribution is observed, 
p(s) ~ s~ a , with a = 2.04 ± 0.04, revealing the absence of a 
characteristic size. The distribution of distances of discon- 
nected loops from the spanning cluster decays exponentially, 
i.e., the loop cloud is mainly localized in the neighborhood of 
the spanning cluster (see Appendix for data). 

To understand how the RGB set emerges, we now consider 
its evolution with the control parameter p. Initially, when the 
fluids R, G, and B start to invade from the boundary, the RGB 
set is empty. As p increases, at a typical value p = p t , two 
fluids for the first time try to invade the same channel in the 
bulk and with increasing p a division surface starts to form 
orthogonal to these channels. If, in addition, any of the four 
edges shared by two neighboring pores of different color is 
also shared by a pore of the third color, an RGB edge emerges. 
Figure 5(b) shows how the total number of RGB nodes, M tot , 
depends on system size at p = p t and at p — p t ± 0.03. While 
above p t results are consistent with the fractal dimension ob- 
served for the final state (p = 1), precisely at p u a negative 
scaling exponent is obtained, M tot ~ L~ l , with t = 0.68 ±0.08. 
This implies that, in the thermodynamic limit, the RGB set is 
empty at p t , while above p t it is fractal of fractal dimension 
dtot- 

If we assume that the formation of RGB nodes at p t is the 
product of two uncorrected processes, namely the formation 
of a dividing surface between two colors, we can show that 
t = d — 2b, where d is the spatial dimension and b is the frac- 
tal dimension of the dividing surfaces at p t (see Appendix). 
For the fractal dimension b of the set connecting two fluids 
at p t , Coniglio has shown in the context of percolation that 
b = 1/ v, where v is the correlation-length critical exponent 
of percolation [8]. In three dimensions, v = 0.8734 ±0.0006 
[9-1 1], giving t = 3 — 2/v ^ 0.71, consistent with our numer- 




FIG. 4. Sketch explaining the presence of disconnected loops in 
the RGB set. We see a horizontal cut through a partitioned medium. 
The two red regions are connected with each other somewhere above 
or below the shown plane. The arrows indicate the existence of a 
loop, as discussed in the main text. On the right, one sees examples 
of cubes of different colors sharing edges. RGB edges are shown as 
thick solid lines. 

ical result. Above the threshold p t the argument leading to the 
expression for t does not hold, since in this regime the inva- 
sion is correlated. Accordingly, the fractal dimension of the 
RGB set is then different from the one of the intersection of 
the two division surfaces. We conjecture that the expression 
for t can be generalized to any dimension d and number of 
different fluids (colors) n, as far as 2 < n < d\ 

t{d,n) = (n-2)d-(n-l)/v. (1) 

For two colors, n = 2, t — — 1/v = —b, as in percolation. In 
contrast for n > 3, given the exact and numerical values for v 
[12], t is positive. Above d = 6, the upper-critical dimension 
of percolation, 1/v = 2, such that t = (n — 2)d -2{n—l). 

We find that M tot scales with the distance to p t as 
^tot ~ (p — Pt)^ T with £7 = 2.0 ±0.3. Therefore, we propose 
the following crossover scaling for the total number of RGB 
nodes: 

M tot (p,L)=L- t G[(p-p t )L e }. (2) 

This scaling behavior of M to t in p and L implies 

= (^tot + O/Cr- I n addition, the scaling function G[x] ful- 
fills G[x] ~ x& for x > 0. The Ansatz in Eq. (2) is confirmed 
by the scaling plot in Fig. 5(c). 



IV. DISCUSSION 

We found a rich scale-free behavior in the partitioning of 
random media through simultaneous invasion by three flu- 
ids. The lines where all three fluids are simultaneously in 
contact form a fractal set, the RGB set, of dimension d tot = 
1.69 ±0.02, while its spanning cluster has dimension d sc — 
1.55 ±0.03. The other clusters are loops and their size fol- 
lows a power-law distribution. At the threshold where two 
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FIG. 5. Scale-free behavior, (a) Number of loop sizes s divided by L dtot as function of s for different lattice sizes L (256 (□), 362 (Q), and 
512 (A)). For intermediate sizes the data follows a power law, p(s) ~ s~ a with a = 2.04 ± 0.04. (b) Size dependence of the total mass M tot 
at and above the threshold pressure p t . At p t (O) the total mass decreases according to a power law, M tot ~ L~\ where t = 0.69 ± 0.08. In 
contrast, above the threshold (□) the slope is given by the fractal dimension d to t- (c) Rescaled total mass M\q\]J as a function of the scaling 
variable (p — pt)L d for different lattice sizes L. Close to the threshold M to t ~ (p — Pt)^ T , where £r = 2.0 ± 0.3. The solid lines are guides to 
the eye. 



fluids first form an interface, the size of the set of RGB nodes 
scales with a negative exponent in the system size. We pro- 
pose a crossover scaling between this exponent and the fractal 
dimension d tot above the threshold. 

For an oil reservoir shared by three companies, our study 
shows how the optimal injection regions scale with the reser- 
voir size and how they are spatially distributed. In the sec- 
ond example, image analysis, our work establishes how the 
number of pixels forming the simultaneous boundary between 
three tissues scales with the image resolution. Besides these 
examples, our results have also implications to other fields. 
Let us consider a chemical reaction that requires three dif- 
ferent reactants each entering a porous medium from another 
side. Supposing that all reactants have the same diffusion 
constant our results identify the disconnected fractal region 
in which the reaction will take place. Finally, knowing the 
properties of the partitioning of a porous medium is also rel- 
evant for planning of waste disposal. Often trash is buried 
under ground, in porous soils, such that its decomposition 
gases spread through pores and fractures [13]. These gases 
will eventually leave the soil and the partitioning of the soil 
determines where this will happen first. In particular, the frac- 
tal RGB set are the disposal regions where the escaping of 
gases occurs, simultaneously, in three regions. 

The exact shape of the RGB set depends on the threshold 
distribution and on the injection areas of the fluids. We also 
tested the partition model with different sets of injection pores, 
namely, (1) division of the six faces of a cubic medium into 
three injection areas, corresponding each to two adjacent faces 
of the cube, (2) injection from three vertical faces of a cube, 
and (3) injection from three edges of the cube, with periodic 
boundary conditions. For all cases, we obtained fractal di- 
mensions consistent with d to t — 1.69 ±0.02. In contrast, for 
the following injection patterns, different values for d tot have 
been obtained: (1) division of the six faces such that each fluid 



is injected from two opposite faces, (2) injection pores uni- 
formly distributed in the cube, with periodic boundary condi- 
tions, and (3) three single injection pores in the cube, also with 
periodic boundaries. These observations suggest that two con- 
ditions on the injection areas, though not necessary, are suffi- 
cient to obtain the RGB fractal dimension reported here. First, 
the injection area of each color must be singly connected. Sec- 
ond, the division of the surface of the medium into these areas 
has to be such that no single fluid can isolate the remaining 
two fluids from each other. 

The reported fractal dimensions were obtained for a uni- 
form and uncorrected distribution of thresholds. It is well- 
known that disorder in soils is typically characterized by spa- 
tial correlations, which can be described by their Hurst ex- 
ponent H. The numerical values of the fractal dimensions 
reported here will in general depend on H [14-17]. Never- 
theless, our qualitative and topological arguments should still 
be applicable. 

Models of discontinuous percolation transitions are a sub- 
ject of recent interest [18-27]. Some of these models lead to 
compact clusters with fractal perimeters [25-27] and in some 
cases with a fractal dimension compatible with the one of di- 
vision surfaces [26]. The simultaneous boundaries between 
three clusters are therefore quite likely related to RGB sets. 

Appendix A: Simulation details 

All numerical results have been obtained with Monte Carlo 
simulations. Random numbers have been generated with the 
algorithm proposed in Ref. [28]. Considering the labeling 
scheme by Newman and Ziff [29, 30], we kept track of the 
color properties as function of the fraction of sampled chan- 
nels p. For Fig. 3, results have been averaged over at least 
2800 realizations. In Fig. 5(a), (b), and (c), results have been 
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averaged over at least 10000, 6400, and 300 realizations, re- 
spectively. Unless indicated otherwise, statistical error bars 
are smaller than the symbol size. The algorithmic procedure 
shares similarities with invasion percolation [31, 32] and the 
fracturing of ranked surfaces [33]. The self- similarity of the 
shortest path in the spanning cluster of the RGB set has been 
confirmed using the yardstick method [34, 35]. 



Appendix B: Scaling exponent t and relation to random 
percolation 

In the context of the scaling behavior at the threshold p u in 
the Section Results, we give the following expression for the 
scaling exponent, 

t = d-2b. (Bl) 

It describes the scaling of the total number of RGB nodes 
with the lattice size as M tot ~ L~ l . This expression can be 
obtained under the assumption that the RGB network, which 
corresponds to the region where the division surfaces merge, 
results (only at p t ) from the intersection of two uncorrected 
dividing surfaces. Suppose that the division surfaces at p t are 
fractals of dimension b, i.e., M sur ~ L b where M sur is the mass 
of the division surface (see explanation below). Then, we can 
make use of a general result for the intersection of two frac- 
tal objects reported in Ref. [36]. Consider the intersection of 
two independent fractals of dimensions df t \ and df^ in spatial 
dimension d, with their centers separated less than the bigger 
radius of gyration of the two. Then the fractal dimension of 
the intersection of both fractals is 

d/,in2 = df,i +df£-d. (B2) 

Therefore, if one considers the intersection of two objects of 
dimension b, one has df^2 — 2b — d. And, since t is defined 
by M tot ~L~\ 

t = -d fAn2 = d-2b. (B3) 

To find expressions for the fractal dimension of the division 
surfaces at the threshold and the value p t of the threshold it- 
self, we apply the relation of our model to random percolation 
[12]. If there would not be the constraint that the fluids cannot 
displace each other, percolation connecting to three sides [37] 
would be recovered. Then, a spanning cluster of invaded pores 
would emerge at a fraction of open channels p — p c -> where p c 
is the bond percolation threshold of the lattice. Therefore, for 
p increasing from zero, the first contact of two fluids in the 
bulk occurs at p = p t = p c , i.e., the contact pressure threshold 
in our model equals the bond percolation threshold of the lat- 
tice. In this work, we have used the bond percolation thresh- 
old of the simple-cubic lattice, p c = 0.2488126 [9, 38]. Let 
us now consider the fractal dimension b of the division sur- 
faces at the threshold. Below p t9 the three fluids are not in 
contact and the set of RGB nodes is empty. Therefore, there 
exists no correlation between the emerging division surfaces 
between the fluids. As a result, the RGB network behaves 



like the intersection between two independent division sur- 
faces. On the other hand, in a non-spanning configuration of 
critical percolation, bonds that, once occupied, would yield 
a spanning cluster are called anti-red bonds and their num- 
ber diverges with the lattice size as L 1 ^ [8, 39], where V is 
the correlation-length critical exponent. Now, opening any of 
the channels contained in the division surfaces of our model 
at p t = p c would give rise to a cluster spanning the entire 
medium. Since, in addition, these surfaces are independent, 
they have the fractal dimension of the anti-red bonds, i.e., 
b = 1/v. For our model, Eq. (B3) yields t = 3 - 2/v « 0.71. 
By applying Eq. (B2), the generalization of the expression for 
t given in Eq. (1) is obtained. Equation (1) applies to our 
model for d = n = 3. There, since the sum of the fractal di- 
mension of two division surfaces is smaller than the spatial 
dimension, the objects are mutually transparent (no intersec- 
tions) and the number of bonds in such surfaces vanishes in 
the thermodynamic limit. Above the threshold since corre- 
lations develop between neighboring surfaces, the consider- 
ations used to derive Eq. (1) for t do not hold. In addition, 
from the numerical values of the fractal dimension, one then 
sees that the RGB network scales differently from the mere 
intersection of uncorrected division surfaces. 



Appendix C: Loops and discretization effects 

In the article, we discuss that, for p = 1, RGB clusters not 
connected to the largest one are loops. An example how they 
emerge is shown in Fig. 4. To understand possible discretiza- 
tion effects on the loops, expected for lattice studies, we ana- 
lyze this example in more detail. Figure 4 shows a sketch of a 
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FIG. 6. CDF of distances of disconnected loops from the spanning 
cluster. The number of loops with distance greater equal k is plotted 
as function of k/L, for different lattice lengths L (64 (□), 128 (OX 
and 256 (A)). 

horizontal cross section through a medium. Since p = 1, the 
medium has been entirely invaded. The dashed circle marks 
the location of an RGB edge belonging to the RGB cluster 
spanning the medium in the direction perpendicular to the pa- 
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per plane, connecting the top and bottom sides of the cube 
in Fig. 2. Inside the dotted circle, one sees the emergence of 
a loop disconnected from the spanning RGB cluster. There, 
the surface dividing the R from the B fluid closely approaches 
the G-R division surface. This leads to the presence of an 
RGB edge at the position indicated by the lower arrow. Fur- 
ther above, where the G and B parts are in contact, the RGB 
network vanishes. Finally, the upper arrow indicates where 
the R-B division surface reappears and another RGB edge 
emerges. Given now the full, three-dimensional medium, not 
only a slice of it but the entire R part is connected and a loop 
emerges in the RGB network. The sketch in Fig. 4 suggests 
a continuum picture, but on the lattice discretization effects 
are possible. Partially developed loops appear as loops with 
a single edge. These effects lead to spurious dangling ends 
attached to the backbone of the largest cluster. Dangling ends 
are not present in the continuum where, in addition to the 



largest cluster, only disconnected loops are expected. 



Appendix D: Distance of disconnected loops from spanning 
cluster 

Figure 6 shows the number of loops separated by a distance 
larger or equal to k from the spanning cluster as a function of 
k/L, where L is the lattice length. The distance of a discon- 
nected loop to the spanning cluster is defined as the minimum 
distance between the loop and the spanning cluster in the dual 
lattice. 
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